"""
CBMPy: Constraint Based Modelling in Python (http://pysces.sourceforge.net/cbm)
============
Copyright (C) 2009-2012 Brett G. Olivier, VU University Amsterdam, Amsterdam, The Netherlands

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>

Author: Brett G. Olivier
Contact email: bgoli@users.sourceforge.net
Last edit: $Author: bgoli $ ($Id: CBModel.py 990 2012-02-15 13:52:11Z bgoli $)

"""

import numpy
from CBDataStruct import StructMatrixLP
from CBVersion import __DEBUG__, __version__

class Fbase(object):
    '''
    Base class for CB Model objects provides *id* and fields
    '''
    id = None
    name = None
    annotation = None
    __TRASH__ = None

    ##  __pre__ = ''
    ##  __post__ = ''

    def getPid(self):
        ##  return self.__pre__ + self.id + self.__post__
        return self.id

    def getName(self):
        return self.name

    def getAnnotations(self):
        return self.annotation

    def setName(self, name):
        self.name = name

    def setAnnotation(self, key, value):
        assert self.annotation != None, '\nThis class has no annotation field'
        self.annotation.update({key : value})

    def setPid(self, id):
        self.id = id


class Model(Fbase):
    '''
    Container for constraint based model, adds methods for manipulating:

     - objectives
     - constraints
     - reactions
     - species
     - N a structmatrix object

    '''
    prefix = ''
    suffix = ''
    objectives = None
    obj_func = None
    activeObjIdx = None
    flux_bounds = None
    reactions = None
    species = None
    N = None
    inputfile_id = ''
    description = ''
    user_constraints = None
    CM = None
    sensitivity = None
    ##  optValue = None
    ##  optType = None
    SOLUTION_STATUS_INT = None
    SOLUTION_STATUS = None
    

    def __init__(self, id):
        '''
        Container for constraint based model.
        
         - *id* a string with the model ID        
        
        Adds methods for manipulating:

         - objectives
         - constraints
         - reactions
         - species
         - N a structmatrix object

        '''        
        self.id = id
        self.objectives = []
        self.flux_bounds = []
        self.reactions = []
        self.species = []
        self.annotation = {}
        self.__TRASH__ = {}

    def addObjective(self, obj, active=False):
        """
        Add an instantiated Objective object to the FBA model
        
        - *obj* an instance of the Objective class
        - *active* [default=False] flag this objective as the active objective (fba.activeObjIdx)
        
        """         
        print 'Adding objective: %s' % obj.id
        self.objectives.append(obj)
        self.obj_func = self.objectives[-1]
        if active:
            self.activeObjIdx = len(self.objectives)-1

    def addFluxBound(self, fluxbound):
        """
        Add an instantiated FluxBound object to the FBA model
        
        - *fluxbound* an instance of the FluxBound class
        
        """                
        RL = self.getFluxBoundByReactionID(fluxbound.getReactionId(), 'lower')
        RU = self.getFluxBoundByReactionID(fluxbound.getReactionId(), 'upper')
        if RL != None and RL.is_bound == fluxbound.is_bound and RL.reaction == fluxbound.reaction:
            print 'Lower bound on reactions %s already exists (adding anyway)' % (fluxbound.getReactionId())
        if RU != None and RU.is_bound == fluxbound.is_bound and RU.reaction == fluxbound.reaction:
            print 'Upper bound on reactions %s already exists (adding anyway)' % (fluxbound.getReactionId())
        if __DEBUG__: print 'Adding fluxbound: %s' % fluxbound.id
        self.flux_bounds.append(fluxbound)

    def addSpecies(self, species):
        """
        Add an instantiated Species object to the FBA model
        
        - *species* an instance of the Species class
        
        """        
        if __DEBUG__: print 'Adding species: %s' % species.id
        self.species.append(species)

    def addReaction(self, reaction):
        """
        Adds a reaction object to the model
        
        - *reaction* an instance of the Reaction class
        
        """
        if __DEBUG__: print 'Adding reaction: %s' % reaction.id
        self.reactions.append(reaction)

    def addUserConstraint(self, id, fluxes=[], operator='=', rhs=0.0):
        """
        Add a user defined constraint to FBA model, this is additional to the automatically determined Stoichiometric constraints.
        
         - *id* user constraint name/id, use `None` for auto-assign
         - *fluxes* a list of (coefficient, reaction id) pairs where coefficient is a float
         - *operator* is one of = > < >= <=
         - *rhs* a float

        """
        if self.user_constraints == None:
            self.user_constraints = {}
        Rids = self.getReactionIds()
        fluxlist = []
        for J in fluxes:
            assert J[1] in Rids, '\n%s is not a valid reaction id' % J[1]
            assert J[1] not in fluxlist, '\nFluxes may only appear once per constraint:\n%s' % fluxes
            fluxlist.append(J[1])
        assert operator in ['>','<','<=','>=','='], '\n%s is not a valid operator' % operator
        rhs = float(rhs)
        if id == None:
            id = 'uConstr%s' % (len(self.user_constraints.keys())+1)
        self.user_constraints.update({id : {'fluxes' : fluxes,
                                              'operator' : operator,
                                              'rhs' : rhs}
                                      })


    def deleteReactionAndBounds(self, reaction_id):
        """
        Delete all reaction and bounds connected to reaction
         
         - *id* a valid reaction id
         
        """
        
        Ridx = None
        Robj = None
        assert reaction_id in self.getReactionIds(), '\nOh Hellooooooooooooo'
        for r in range(len(self.reactions)-1, -1, -1):
            if self.reactions[r].getPid() == reaction_id:
                Ridx = reaction_id
                Robj = self.reactions.pop(r)
        Bounds = []
        for b in range(len(self.flux_bounds)-1, -1, -1):
            if self.flux_bounds[b].reaction == reaction_id:
                Bounds.append(self.flux_bounds.pop(b))
        print 'Deleting reaction %s and %s associated bounds' % (Ridx, len(Bounds))
        self.__TRASH__.update({Ridx : {'react' : Robj, 'bnds' : Bounds}})

    def deleteObjective(self, objective_id):
        """
        Delete objective function:
        
         *objective_id* the id of the objective function. If objective_id is given  as 'active' then the active objective is deleted.
        
        """
        if objective_id == 'active':
            objective_id = self.getActiveObjective().getPid()
        
        assert objective_id in self.getObjectiveIds(), '\nNo ....'
        for o in range(len(self.objectives)-1, -1, -1):
            if self.objectives[o].getPid() == objective_id:
                Oidx = objective_id
                Oobj = self.objectives.pop(o)
        print 'Deleting objective %s' % (Oidx)
        self.__TRASH__.update({Oidx : Oobj})
        
    def undeleteObjective(self, objective_id):
        """
        Undeltes a deleted objective function:
        
         - *objective_id* the id of an objeective function
         
        """
        
        assert self.__TRASH__.has_key(objective_id), '\nNo deleted object of with this id'
        self.addObjective(self.__TRASH__[objective_id])

        
    def deleteBoundsForReactionId(self, reaction_id):
        """
        Delete all bounds connected to reaction, reaction_id
         
         - *reaction_id* a valid reaction id
         
        """      
        ##  Ridx = None
        ##  Robj = None
        Bounds = []
        for b in range(len(self.flux_bounds)-1, -1, -1):
            if self.flux_bounds[b].reaction == reaction_id:
                Bounds.append(self.flux_bounds.pop(b))
        print 'Deleting %s bounds associated with reaction %s' % (len(Bounds), reaction_id)
        ##  self.__TRASH__.update({Ridx : {'react' : Robj, 'bnds' : Bounds}})


    def undeleteReactionAndBounds(self, reaction_id):
        """
        Undelete a reaction and bounds deleted with the **deleteReactionAndBounds** method
        
         - *id* a deleted reaction id
         
        Please note this method is still experimental ;-)
        
        """
        assert self.__TRASH__.has_key(reaction_id), '\nOops I did it again ...'
        Ridx = reaction_id
        self.addReaction(self.__TRASH__[reaction_id]['react'])
        for b in self.__TRASH__[reaction_id]['bnds']:
            self.addFluxBound(b)
        self.__TRASH__.pop(reaction_id)
        
        print 'Undeleting reaction: %s' % reaction_id

    def changeAllFluxBoundsWithValue(self, old, new):
        """
        Replaces all flux bounds with value "old" with a new value "new":

         - *old* value
         - *new* value

        """
        cntr = 0
        for b in range(len(self.flux_bounds)-1,-1,-1):
            if self.flux_bounds[b].value == old:
                ##  db = fba.flux_bounds.pop(b)
                self.flux_bounds[b].value = new
                cntr += 1
        print '\nReplaced %s \"%s\" bounds with \"%s\"' % (cntr, old, new)

    def deleteAllFluxBoundsWithValue(self, value):
        """
        Delete all flux bounds which have a specified value:

         - *value* the value of the flux bound(s) to delete

        """
        cntr = 0
        for b in range(len(self.flux_bounds)-1,-1,-1):
            if self.flux_bounds[b].value == value:
                db = self.flux_bounds.pop(b)
                cntr += 1
        print '\nDeleted %s \"%s\" bounds' % (cntr, value)

# TODO: undelete Species
    def deleteSpecies(self, sid):
        """
        Deletes a species object with id
        
        - *sid* the species id
        
        """
        s = self.species.pop(self.getSpeciesIds().index(sid))

    def deleteNonReactingSpecies(self, simulate=True):
        """
        Deletes all species that are not reagents (do not to take part in a reaction). 
        *Warning* this deletion is permanent and greedy (not selective). Returns a list of (would be) deleted species
        
         - *simulate* [default=True] only return a list of the speciesId's that would have been deleted if False
        
        """
        active_reagents = []
        for R in self.reactions:
            for S in R.getReagentRefs():
                if S not in active_reagents:
                    active_reagents.append(S)
            if __DEBUG__: print active_reagents
        if not simulate:
            print 'Deleting non-reactive species',
        deleted_species = []
        for S in range(len(self.species)-1,-1,-1):
            if self.species[S].getPid() not in active_reagents:
                deleted_species.append(self.species[S].getPid())
                if simulate:
                    pass
                else:
                    self.deleteSpecies(self.species[S].getPid())
                if __DEBUG__: print self.getSpeciesIds()
                print '.',
        print ' '
        return deleted_species

    def getReaction(self, rid):
        """
        Returns a reaction object with *id*
        
         - *rid* reaction ID
        
        """
        for r in self.reactions:
            if r.getPid() == rid:
                return r

    def getSpecies(self, sid):
        """
        Returns a species object with *sid*
        
         - *sid* a specied ID
        
        """
        for s in self.species:
            if s.getPid() == sid:
                return s

    def getReactionBounds(self, rid):
        """
        Get the bounds of a reaction, returns a tuple of rid, lowerbound value, upperbound value and equality value (None means bound does not exist).
        
         - *rid* the reaction ID
        
        """
        lb = ub = eq = None
        lb = self.getFluxBoundByReactionID(rid, 'lower')
        ub = self.getFluxBoundByReactionID(rid, 'upper')
        eq = self.getFluxBoundByReactionID(rid, 'equality')
        if lb != None:
            if numpy.isinf(lb.value) or numpy.isreal(lb.value):
                lb = lb.value
            else:
                lb = float(lb)
        if ub != None:
            if numpy.isinf(ub.value) or numpy.isreal(ub.value):
                ub = ub.value
            else:
                ub = float(ub)
        if eq != None:
            if numpy.isinf(eq.value) or numpy.isreal(eq.value):
                eq = eq.value
            else:
                eq = float(eq)
        return rid, lb, ub, eq

    def getReactionLowerBound(self, rid):
        """
        Returns the lower bound of a reaction (it it exists) or None
        
         - *rid* the reaction ID
        
        """
        lb = eq = None
        lb = self.getFluxBoundByReactionID(rid, 'lower')
        if lb != None:
            if numpy.isinf(lb.value) or numpy.isreal(lb.value):
                lb = lb.value
            else:
                lb = float(lb)
        else:
            eq = self.getFluxBoundByReactionID(rid, 'equality')
            if eq != None:
                print '\nLower bound defined as an equality'
                if numpy.isinf(eq.value) or numpy.isreal(eq.value):
                    lb = eq.value
                else:
                    lb = float(eq)
        return lb

    def getReactionUpperBound(self, rid):
        """
        Returns the upper bound of a reaction (it it exists) or None
        
         - *rid* the reaction ID
                
        """
        ub = eq = None
        ub = self.getFluxBoundByReactionID(rid, 'upper')
        if ub != None:
            if numpy.isinf(ub.value) or numpy.isreal(ub.value):
                ub = ub.value
            else:
                ub = float(ub)
        else:
            eq = self.getFluxBoundByReactionID(rid, 'equality')
            if eq != None:
                print '\nUpper bound defined as an equality'
                if numpy.isinf(eq.value) or numpy.isreal(eq.value):
                    ub = eq.value
                else:
                    ub = float(eq)
        return ub

    ##  def getBoundByName(self, rid, bound):
        ##  """
        ##  Return a FluxBound instance. Note this is an old name for the newer preferred method: `getFluxBoundByReactionID`
        
         ##  - *rid* the reaction ID
         ##  - *bound* the bound: 'upper','lower','equal'
                
        ##  """
        ##  print 'Deprecated: use *getFluxBoundByReactionID*'
        ##  return self.getFluxBoundByReactionID(rid, bound)

    def getFluxBoundByReactionID(self, rid, bound):
        """
        Returns a FluxBound instance
        
         - *rid* the reaction ID
         - *bound* the bound: 'upper','lower','equal'
                
        """
        c2 = [c for c in self.flux_bounds if c.reaction == rid and c.is_bound == bound]
        if len(c2) >= 1:
            if __DEBUG__: print c2[0].id, c2[0].is_bound, c2[0].reaction, c2[0].operation, c2[0].value
            return c2[0]
        else:
            return None

    def getSpeciesIds(self, hasstring=None):
        if hasstring == None:
            return [s.getPid() for s in self.species]
        else:
            return [s.getPid() for s in self.species if hasstring in s.getPid()]

    def getReactionIds(self, hasstring=None):
        if hasstring == None:
            return [r.getPid() for r in self.reactions]
        else:
            return [r.getPid() for r in self.reactions if hasstring in r.getPid()]

    def getReactionNames(self, hasstring=None):
        if hasstring == None:
            return [r.name for r in self.reactions]
        else:
            return [r.name for r in self.reactions if hasstring in r.name]

    def getFluxBoundIds(self):
        return [s.getPid() for s in self.flux_bounds]

    def getObjectiveIds(self):
        return [s.getPid() for s in self.objectives]

    def getActiveObjective(self):
        """
        Returns the active objective object.
        
        """
        return self.objectives[self.activeObjIdx]
        

    def setActiveObjective(self, objId):
        idz = [o.id for o in self.objectives]
        if objId in idz:
            self.activeObjIdx = idz.index(objId)
            self.obj_func = self.objectives[self.activeObjIdx]
            print 'Objective %s (idx=%s) set active' % (objId, self.activeObjIdx)
        else:
            print 'Invalid ObjectiveId: %s' % objId

    def setBoundValueByName(self, r_id, value, bound):
        """
        Deprecated use setReactionBound
        
        Set a reaction bound
         - *r_id* the reactions id
         - *value* the new value
         - *bound* this is either 'lower' or 'upper'

        """
        print "Deprecated: use *setReactionBound*"
        self.setReactionBound(r_id, value, bound)
        

    def setReactionBound(self, r_id, value, bound):
        """
        Set a reaction bound

         - *r_id* the reactions id
         - *value* the new value
         - *bound* this is either 'lower' or 'upper'

        """
        c2 = self.getFluxBoundByReactionID(r_id, bound)
        # changed to no str() casting
        if c2 != None:
            if numpy.isreal(value):
                c2.value = value
            elif type(value) == unicode or type(value) == str:
                c2.value = float(value)            
            elif numpy.isinf(value):
                c2.value = value
            else:
                c2.value = float(value)            
            if __DEBUG__: print c2.reaction, c2.operation, c2.value
        else:
            print '%s is not a valid reaction name or \'%s\' bound does not exist ' % (r_id,bound)
        
    def setReactionLowerBound(self, r_id, value):
        """
        Set a reactions lower bound (if it exists)

         - *r_id* the reactions id
         - *value* the new value

        """        
        self.setReactionBound(r_id, value, 'lower')

    def setReactionUpperBound(self, r_id, value):
        """
        Set a reactions upper bound (if it exists)

         - *r_id* the reaction id
         - *value* the new value

        """
        self.setReactionBound(r_id, value, 'upper')

    def setPrefix(self, prefix, target):
        """
        This is alpha stuff, target can be:

         - 'species'
         - 'reactions'
         - 'constraints'
         - 'objectives'
         - 'all'

        """
        self.prefix = prefix
        ##  self.id = prefix+self.id
        if target == 'species' or target == 'all':
            for s in self.species:
                s.id = prefix+s.id
            if hasattr(self,'N'):
                if hasattr(self.N,'row'):
                    for s in range(len(self.N.row)):
                        self.N.row[s] = prefix+self.N.row[s]
        if target == 'reactions' or target == 'all':
            for s in self.reactions:
                s.id = prefix+s.id
                for r in s.reagents:
                    r.id = prefix+r.id
                    r.species_ref = prefix+r.species_ref
            if hasattr(self,'N'):
                if hasattr(self.N,'col'):
                    for r in range(len(self.N.col)):
                        self.N.col[r] = prefix+self.N.col[r]
        if target == 'constraints' or target == 'all':
            for s in self.flux_bounds:
                s.id = prefix+s.id
                s.reaction = prefix+s.reaction
        if target == 'objectives' or target == 'all':
            for s in self.objectives:
                s.id = prefix+s.id
                for f in s.fluxObjectives:
                    f.id = prefix+f.id
                    f.reaction = prefix+f.reaction

    def setSuffix(self, suffix, target):
        """
        This is alpha stuff, target can be:

         - 'species'
         - 'reactions'
         - 'constraints'
         - 'objectives'
         - 'all'

        """
        self.suffix = suffix
        ##  self.id = self.id+suffix
        if target == 'species' or target == 'all':
            for s in self.species:
                s.id = s.id+suffix
            if hasattr(self,'N'):
                if hasattr(self.N,'row'):
                    for s in range(len(self.N.row)):
                        self.N.row[s] = self.N.row[s]+suffix
        if target == 'reactions' or target == 'all':
            for s in self.reactions:
                s.id = s.id+suffix
                for r in s.reagents:
                    r.id = r.id+suffix
                    r.species_ref = r.species_ref+suffix
            if hasattr(self,'N'):
                if hasattr(self.N,'col'):
                    for r in range(len(self.N.col)):
                        self.N.col[r] = self.N.col[r]+suffix
        if target == 'constraints' or target == 'all':
            for s in self.flux_bounds:
                s.id = s.id
                s.reaction = s.reaction+suffix
        if target == 'objectives' or target == 'all':
            for s in self.objectives:
                s.id = s.id+suffix
                for f in s.fluxObjectives:
                    f.id = f.id+suffix
                    f.reaction = f.reaction+suffix
    

    def setObjectiveFlux(self, reaction_id, coefficient=1, osense='maximize', delete_objflx=True):
        """
        Set single target reaction flux for the current active objective function.

         - *reaction_id* a string containing a reaction id
         - *coefficient* [default=1] an objective flux coefficient 
         - *osense* the optimization sense must be **maximize** or **minimize**
         - *delete_objflx* [default=True] delete all existing fluxObjectives in the active objective function
        """
        assert reaction_id in self.getReactionIds(), '\n%s is not a reaction\n%s' % (reaction_id, self.getReactionIds())
        osense = osense.lower()
        if osense == 'max': osense = 'maximize'
        if osense == 'min': osense = 'minimize'        
        if osense in ['maximise', 'minimise']:
            osense = osense.replace('se','ze')
        assert osense in ['maximize', 'minimize'], "\nosense must be ['maximize', 'minimize'] not %s" % osense
        if delete_objflx:
            self.objectives[self.activeObjIdx].fluxObjectives = []        
        FO = FluxObjective(reaction_id+'_fluxobj', reaction_id, coefficient)
        self.objectives[self.activeObjIdx].addFluxObjective(FO)
        self.objectives[self.activeObjIdx].operation = osense

    def sortReactionsById(self):
        """
        Sorts the reactions by Reaction.id uses the python string sort
        
        """
        
        reactions2 = []
        reactId = self.getReactionIds()
        reactId.sort()
        for r in reactId:
            reactions2.append(self.getReaction(r))
        self.reactions = reactions2
        del reactions2

    def sortSpeciesById(self):
        """
        Sorts the reactions by Reaction.id uses the python string sort
        
        """
        raise NotImplementedError, "\nCall developer about this"

    def getExchangeReactions(self):
        """
        Returns reaction instances where the 'is_exchange' attribute set to True. This is by default
        reactions that contain a boundary species.
        """
        return [r for r in self.reactions if r.is_exchange]

    def getExchangeReactionIds(self):
        """
        Returns id's of reactions where the 'is_exchange' attribute set to True. This is by default
        reactions that contain a boundary species.
        """
        return [r.getPid() for r in self.reactions if r.is_exchange]

    def getReactionValues(self, only_exchange=False):
        """
        Returns a dictionary of ReactionID : ReactionValue pairs:
        
         - *only_exchange* [default=False] only return the reactions labelled as exchange
         
        """
        out = {}
        reacts = []
        if only_exchange:
            reacts = self.getExchangeReactionIds()
        else:
            reacts = self.getReactionIds()
        for rid in reacts:
            out.update({rid : self.getReaction(rid).value})
        return out

    def buildStoichMatrix(self):
        """
        Build the stoichiometric matrix N and additional constraint matrix CN (if required)
        """
        if __DEBUG__: print 'Species:', self.getSpeciesIds()
        var_spec = [s for s in self.species if not s.is_boundary]
        var_spec_id = [s.getPid() for s in self.species if not s.is_boundary]
        reac_id = self.getReactionIds()
        num_col = len(self.reactions)
        num_row = len(var_spec)
        if __DEBUG__: print 'N-dimension = (%s, %s)' % (num_row, num_col)
        N = numpy.zeros((num_row, num_col))
        for c in range(num_col):
            if __DEBUG__: print self.reactions[c].getPid()
            if __DEBUG__: print self.reactions[c].getStoichiometry()
            for reag in self.reactions[c].getStoichiometry():
                if reag[1] in var_spec_id:
                    if __DEBUG__: print '%s: setting reagent %s to %s (idx=%s,%s)' % (self.reactions[c].getPid(), reag[1], reag[0], var_spec_id.index(reag[1]), c)
                    N[var_spec_id.index(reag[1]), c] = reag[0]
        N = StructMatrixLP(N, range(num_row), range(num_col), row=var_spec_id, col=reac_id)
        self.N = N

        # build and append additional constraint matric
        if self.user_constraints != None:
            crows = self.user_constraints.keys()
            crows.sort()
            ccols = reac_id
            cnum_col = len(ccols)
            cnum_row = len(crows)
            CM = numpy.zeros((cnum_row,cnum_col))
            CRHS = numpy.zeros(cnum_row)
            Coperators = ['=']*cnum_row
            for cs in range(cnum_row):
                for flx in self.user_constraints[crows[cs]]['fluxes']:
                    CM[cs,ccols.index(flx[1])] = float(flx[0])
                    Coperators[cs] = self.user_constraints[crows[cs]]['operator']
                    CRHS[cs] = self.user_constraints[crows[cs]]['rhs']
            self.CM = StructMatrixLP(CM, range(cnum_row), range(cnum_col), row=crows, col=ccols,\
                            rhs=CRHS, operators=Coperators)





class Objective(Fbase):
    '''An objective function'''
    fluxObjectives = None
    ##  fluxObjectiveNames = None
    operation = None
    value = None
    solution = None

    def __init__(self, id, operation):
        self.id = id
        self.operation = operation
        self.fluxObjectives = []

    def addFluxObjective(self, fobj):
        self.fluxObjectives.append(fobj)

    def createFluxObjectives(self, fluxlist=[]):
        """
        Create and add flux objective objects to this objective function.
        
         - *fluxlist* a list of one or more ('coefficient','Reaction_id') pairs
        
        """
        for J in fluxlist:
            if J[1] not in self.getFluxObjectiveReactions():
                self.addFluxObjective(FluxObjective(J[1]+'_fobj', J[1], J[0]))
            else:
                print 'Objective %s already contains flux %s ... skipping!' % (self.getPid(), J[1])

    def getFluxObjectiveIDs(self):
        return [f.getPid() for f in self.fluxObjectives]

    def getFluxObjectiveReactions(self):
        return [f.reaction for f in self.fluxObjectives]

    def getValue(self):
        """
        Returns the current value of the attribute (input/solution)
        """
        return self.value

    def setValue(self, value):
        """
        Sets the attribute ''value''
        """
        self.value = value


class FluxObjective(Fbase):
    '''
    A weighted flux that appears in an objective function

    NOTE: reaction is a string containing a reaction id
    '''
    reaction = None
    coefficient = None

    def __init__(self, id, reaction, coefficient=1):
        self.id = id
        self.reaction = reaction
        self.coefficient = coefficient
        self.annotation = {}

    def getReactionId(self):
        return self.reaction
        
    def getCoefficient(self):
        return self.coefficient

class FluxBound(Fbase):
    '''A reaction fluxbound'''
    reaction = None
    operation = None
    value = None
    is_bound = None

    def __init__(self, id, reaction, operation, value):
        self.id = id
        self.reaction = reaction
        self.operation = operation
        self.value = value
        if operation in ['greater','greaterEqual']:
            self.is_bound = 'lower'
        elif operation in ['less','lessEqual']:
            self.is_bound = 'upper'
        else:
            self.is_bound = 'equality'
        self.annotation = {}

    def getReactionId(self):
        return self.reaction

    def getValue(self):
        """
        Returns the current value of the attribute (input/solution)
        """
        return self.value

    def setValue(self, value):
        """
        Sets the attribute ''value''
        """
        self.value = value


class Reaction(Fbase):
    '''Holds reaction information'''
    reagents = None
    reversible = None
    is_exchange = False
    value = None
    reduced_cost = None
    compartment = None

    def __init__(self, id, name=None, reversible=True):
        self.id = id
        self.name = name
        self.reagents = []
        self.reversible = reversible
        self.annotation = {}
        self.__TRASH__ = {}

    def addReagent(self, reag):
        self.reagents.append(reag)

    def getReagentIds(self):
        return [r.getPid() for r in self.reagents]

    def getReagentRefs(self):
        return [r.species_ref for r in self.reagents]

    def getReagent(self, species):
        print 'Please use R.getReagentWithSpeciesRef(species) instead!\n'
        rgnt = [r for r in self.reagents if r.species_ref == species]
        if len(rgnt) == 0:
            return None
        else:
            return rgnt[0]

    def getValue(self):
        """
        Returns the current value of the attribute (input/solution)
        """
        return self.value

    def setValue(self, value):
        """
        Sets the attribute ''value''
        """
        self.value = value
        
    def getReagentWithSpeciesRef(self, species):
        rgnt = [r for r in self.reagents if r.species_ref == species]
        if len(rgnt) == 0:
            return None
        else:
            return rgnt[0]

    def getStoichiometry(self):
        return [(r.coefficient, r.species_ref) for r in self.reagents]

    def deleteReagent(self, species):
        rids = self.getReagentRefs()
        assert species in rids, '\nWooooaaa boy wooaaa'
        self.__TRASH__.update({self.getPid()+species : self.reagents.pop(rids.index[species])})

    def undeleteReagent(self, species):
        assert self.__TRASH__.has_key(species), '\nAha, yes, sure, maybe, perhaps ...'
        self.reagents.append(self.__TRASH__.pop(self.getPid()+species))


class Species(Fbase):
    '''Holds species information'''
    chemFormula = None
    charge = None
    value = None
    is_boundary = False
    compartment = None
    reagent_of = None

    def __init__(self, id, boundary=False, name=None, value=None,\
                    compartment=None, charge=None, chemFormula=None):
        self.id = id
        self.name = name
        self.value = value
        self.is_boundary = boundary
        self.compartment = compartment
        self.charge = None
        self.chemFormula = None
        self.reagent_of = []
        self.annotation = {}

    def getValue(self):
        """
        Returns the current value of the attribute (input/solution)
        """
        return self.value

    def setValue(self, value):
        """
        Sets the attribute ''value''
        """
        self.value = value
        

class Reagent(Fbase):
    '''
    Has a reactive species id and stoichiometric coefficient:
     - negative = substrate
     - positive = product
     - species_ref = id = a reference to a species obj
    '''
    coefficient = None
    role = None
    species_ref = None

    def __init__(self, id, species_ref, coef):
        self.id = id
        self.species_ref = species_ref
        self.coefficient = coef
        if coef > 0.0:
            self.role = 'product'
        elif coef < 0.0:
            self.role = 'substrate'
        else:
            raise RuntimeError, '\nZero coefficients are NOT cool: {%s}%s!' % (coef, id)
        self.annotation = {}

if __name__ == '__main__':
    s1 = Species('S0', name='S 0', boundary=False, value=0.0,\
        compartment='cytosol', charge=-2, chemFormula='C2H4O')
    s2 = Species('S1', name='S 1', boundary=False, value=0.0,\
        compartment='cytosol', charge=0, chemFormula='C2H4O')
    gs = Reagent('rs1', -1)
    gp = Reagent('rp1', 1)

    R = Reaction('R1',reversible=True)
    R.addReagent(gs)
    R.addReagent(gp)

    c1 = FluxBound('r1c1', R, '>=', -10)
    c2 = FluxBound('r1c2', R, '<=', 10)

    FJ1 = Reaction('FluxObj1', reversible=True)
    FJ2 = Reaction('FluxObj2', reversible=True)
    FJ3 = Reaction('FluxObj3', reversible=True)

    foj1 = FluxObjective('foj1', FJ1, 1)
    foj2 = FluxObjective('foj2', FJ2, 2)
    foj3 = FluxObjective('foj3', FJ3, 1)

    O1 = Objective('objective1', 'maximize')
    O1.addFluxObjective(foj3)

    O2 = Objective('objective2', 'maximize')
    O2.addFluxObjective(foj1)
    O2.addFluxObjective(foj2)

    print O1.getFluxObjectiveIDs()
    print O2.getFluxObjectiveIDs()

    mod = Model('TestModel')
    mod.addObjective(O1, active=False)
    print mod.activeObjIdx
    mod.addObjective(O2, active=True)
    print mod.activeObjIdx
    mod.setActiveObjective('objective1')
    print mod.activeObjIdx

    mod.addFluxBound(c1)
    mod.addFluxBound(c2)
    mod.addSpecies(s1)
    mod.addSpecies(s2)
    mod.addReaction(R)
    mod.addReaction(FJ1)
    mod.addReaction(FJ2)
    mod.addReaction(FJ3)

    print mod.species
    print mod.getSpeciesIds()
    print mod.reactions
    print mod.getReactionIds
    print mod.flux_bounds
    print mod.getFluxBoundIds()
    print mod.objectives
    print mod.getObjectiveIds()
    print R.getReagentRefs()
    print R.getStoichiometry()
    print c1.getReactionId()
